#include "settings.h"
#include "cylinder_tct.h"
#include "fem/trianglefem.h"
#include "fem/triangulation.h"
#include "geometry/2d/vector.h"
#include "geometry/2d/geometry.h"
#include "geometry/2d/lineparam.h"
#include "approximation/interpolation.h"

CylinderTCt::CylinderTCt() : SolverTCt()
{
    //****************************************************************************************************************************
    // Инициализация поисковика по сетке.
    //****************************************************************************************************************************
    m_NetSearcher = new fem::TriangleNetSearch(30, 60, 0.0025, 0.0025);

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    QVariantList r = Settings::instance()->setting("Crucible/R").toList();
    QVariantList z = Settings::instance()->setting("Crucible/Z").toList();
    QVariantList n = Settings::instance()->setting("Crucible/N").toList();

    for(int i = 0; i < r.size(); ++i)
    {
        m_R[i] = r[i].toDouble();
        m_Z[i] = z[i].toDouble();
    }

    for(int i = 0; i < n.size(); ++i)
        m_N[i] = n[i].toDouble();

    m_L[0] = geometry::_2d::distance(geometry::_2d::Point(m_R[0], m_Z[0]), geometry::_2d::Point(m_R[1], m_Z[1]));
    m_L[1] = geometry::_2d::distance(geometry::_2d::Point(m_R[0], m_Z[0]), geometry::_2d::Point(m_R[2], m_Z[2]));
    m_L[2] = geometry::_2d::distance(geometry::_2d::Point(m_R[1], m_Z[1]), geometry::_2d::Point(m_R[3], m_Z[3]));
    m_L[3] = geometry::_2d::distance(geometry::_2d::Point(m_R[2], m_Z[2]), geometry::_2d::Point(m_R[3], m_Z[3]));
    m_L[4] = geometry::_2d::distance(geometry::_2d::Point(m_R[2], m_Z[2]), geometry::_2d::Point(m_R[4], m_Z[4]));
    m_L[5] = geometry::_2d::distance(geometry::_2d::Point(m_R[3], m_Z[3]), geometry::_2d::Point(m_R[5], m_Z[5]));
    m_L[6] = geometry::_2d::distance(geometry::_2d::Point(m_R[1], m_Z[1]), geometry::_2d::Point(m_R[6], m_Z[6]));
    m_L[7] = geometry::_2d::distance(geometry::_2d::Point(m_R[4], m_Z[4]), geometry::_2d::Point(m_R[5], m_Z[5]));
    m_L[8] = geometry::_2d::distance(geometry::_2d::Point(m_R[5], m_Z[5]), geometry::_2d::Point(m_R[6], m_Z[6]));
    m_L[9] = geometry::_2d::distance(geometry::_2d::Point(m_R[4], m_Z[4]), geometry::_2d::Point(m_R[7], m_Z[7]));
    m_L[10] = geometry::_2d::distance(geometry::_2d::Point(m_R[5], m_Z[5]), geometry::_2d::Point(m_R[8], m_Z[8]));
    m_L[11] = geometry::_2d::distance(geometry::_2d::Point(m_R[6], m_Z[6]), geometry::_2d::Point(m_R[10], m_Z[10]));
    m_L[12] = geometry::_2d::distance(geometry::_2d::Point(m_R[9], m_Z[9]), geometry::_2d::Point(m_R[8], m_Z[8]));
    m_L[13] = geometry::_2d::distance(geometry::_2d::Point(m_R[8], m_Z[8]), geometry::_2d::Point(m_R[10], m_Z[10]));
    m_L[14] = geometry::_2d::distance(geometry::_2d::Point(m_R[7], m_Z[7]), geometry::_2d::Point(m_R[11], m_Z[11]));
    m_L[15] = geometry::_2d::distance(geometry::_2d::Point(m_R[7], m_Z[7]), geometry::_2d::Point(m_R[9], m_Z[9]));
    m_L[16] = geometry::_2d::distance(geometry::_2d::Point(m_R[9], m_Z[9]), geometry::_2d::Point(m_R[10], m_Z[10]));
    m_L[17] = geometry::_2d::distance(geometry::_2d::Point(m_R[9], m_Z[9]), geometry::_2d::Point(m_R[12], m_Z[12]));
    m_L[18] = geometry::_2d::distance(geometry::_2d::Point(m_R[11], m_Z[11]), geometry::_2d::Point(m_R[12], m_Z[12]));

    //****************************************************************************************************************************
    // Создание конечноэлементной сетки.
    //****************************************************************************************************************************
    fem::Triangulation triangulation(QStringList() << "T");

    int basis_points[13];
    basis_points[0] =  triangulation.addBasisPoint(geometry::_2d::Point(m_R[0],m_Z[0]));
    basis_points[1] =  triangulation.addBasisPoint(geometry::_2d::Point(m_R[1],m_Z[1]));
    basis_points[2] =  triangulation.addBasisPoint(geometry::_2d::Point(m_R[2],m_Z[2]));
    basis_points[3] =  triangulation.addBasisPoint(geometry::_2d::Point(m_R[3],m_Z[3]));
    basis_points[4] =  triangulation.addBasisPoint(geometry::_2d::Point(m_R[4],m_Z[4]));
    basis_points[5] =  triangulation.addBasisPoint(geometry::_2d::Point(m_R[5],m_Z[5]));
    basis_points[6] =  triangulation.addBasisPoint(geometry::_2d::Point(m_R[6],m_Z[6]));
    basis_points[7] =  triangulation.addBasisPoint(geometry::_2d::Point(m_R[7],m_Z[7]));
    basis_points[8] =  triangulation.addBasisPoint(geometry::_2d::Point(m_R[8],m_Z[8]));
    basis_points[9] =  triangulation.addBasisPoint(geometry::_2d::Point(m_R[9],m_Z[9]));
    basis_points[10] = triangulation.addBasisPoint(geometry::_2d::Point(m_R[10],m_Z[10]));
    basis_points[11] = triangulation.addBasisPoint(geometry::_2d::Point(m_R[11],m_Z[11]));
    basis_points[12] = triangulation.addBasisPoint(geometry::_2d::Point(m_R[12],m_Z[12]));

    int border_segments[19];
    border_segments[0] =  triangulation.addBorderSegment(basis_points[0], basis_points[1], m_N[0]);
    border_segments[1] =  triangulation.addBorderSegment(basis_points[0], basis_points[2], m_N[1]);
    border_segments[2] =  triangulation.addBorderSegment(basis_points[1], basis_points[3], m_N[2]);
    border_segments[3] =  triangulation.addBorderSegment(basis_points[2], basis_points[3], m_N[3]);
    border_segments[4] =  triangulation.addBorderSegment(basis_points[2], basis_points[4], m_N[4]);
    border_segments[5] =  triangulation.addBorderSegment(basis_points[3], basis_points[5], m_N[5]);
    border_segments[6] =  triangulation.addBorderSegment(basis_points[1], basis_points[6], m_N[6]);
    border_segments[7] =  triangulation.addBorderSegment(basis_points[4], basis_points[5], m_N[7]);
    border_segments[8] =  triangulation.addBorderSegment(basis_points[5], basis_points[6], m_N[8]);
    border_segments[9] =  triangulation.addBorderSegment(basis_points[4], basis_points[7], m_N[9]);
    border_segments[10] = triangulation.addBorderSegment(basis_points[5], basis_points[8], m_N[10]);
    border_segments[11] = triangulation.addBorderSegment(basis_points[6], basis_points[10],m_N[11]);
    border_segments[12] = triangulation.addBorderSegment(basis_points[7], basis_points[8], m_N[12]);
    border_segments[13] = triangulation.addBorderSegment(basis_points[8], basis_points[10],m_N[13]);
    border_segments[14] = triangulation.addBorderSegment(basis_points[7], basis_points[11],m_N[14]);
    border_segments[15] = triangulation.addBorderSegment(basis_points[7], basis_points[9], m_N[15]);
    border_segments[16] = triangulation.addBorderSegment(basis_points[9], basis_points[10],m_N[16]);
    border_segments[17] = triangulation.addBorderSegment(basis_points[9], basis_points[12],m_N[17]);
    border_segments[18] = triangulation.addBorderSegment(basis_points[11],basis_points[12],m_N[18]);

    triangulation.addBorderSegmentBCInfo("T", border_segments[0], 1, 0);
    triangulation.addBorderSegmentBCInfo("T", border_segments[6], 1, 1);
    triangulation.addBorderSegmentBCInfo("T", border_segments[11], 1, 1);
    triangulation.addBorderSegmentBCInfo("T", border_segments[18], 1, m_TRod != 0 ? 2 : 1);

    triangulation.addQuadrangleSubdomain(border_segments[3],  border_segments[2],  border_segments[0],  border_segments[1], 0);
    triangulation.addQuadrangleSubdomain(border_segments[7],  border_segments[5],  border_segments[3],  border_segments[4], 1);
    triangulation.addQuadrangleSubdomain(border_segments[8],  border_segments[6],  border_segments[2],  border_segments[5], 0);
    triangulation.addQuadrangleSubdomain(border_segments[12], border_segments[10], border_segments[7],  border_segments[9], 3);
    triangulation.addQuadrangleSubdomain(border_segments[13], border_segments[11], border_segments[8],  border_segments[10],0);
    triangulation.addQuadrangleSubdomain(border_segments[16], border_segments[13], border_segments[12], border_segments[15],0);
    triangulation.addQuadrangleSubdomain(border_segments[18], border_segments[17], border_segments[15], border_segments[14],0);

    m_State = SolverTCt::NotStarted;

    m_Net = triangulation.triangulate();
    m_NetSearcher->analyzeNet(m_Net);

    //****************************************************************************************************************************
    //Формирование информации о краевых условиях.
    //****************************************************************************************************************************
    triangulation.getDirichletBoundaryInfo("T", m_TopTBC, m_TypeTBC);

    //****************************************************************************************************************************
    // Формирование начальных значений для векторов решения.
    //****************************************************************************************************************************
    m_T = linal::Vector(m_Net.numOfPoints(), m_TMax);
    m_TOldTime = linal::Vector(m_Net.numOfPoints(), m_TMax);

    //****************************************************************************************************************************
    // Формирование портрета СЛАУ. Выделение памяти под СЛАУ.
    //****************************************************************************************************************************
    m_F = linal::Vector(m_Net.numOfPoints());

    QVector<int> igl, jgl, igu, nju;
    fem::triangles::formMatrixProfile(m_Net, igl, jgl, igu, nju);
    m_A = linal::CSRMatrix(m_Net.numOfPoints(), igl, jgl, igu, nju);

    foreach(linal::SLAESolver* solver, m_Solvers)
        solver->adaptForSLAE(m_A.dimension());

    //****************************************************************************************************************************
    // Графит.
    //****************************************************************************************************************************
    m_Materials.append(Material(Material::Solid));
    m_Materials[0].setProperty(Material::Diffusivity, 0);
    m_Materials[0].setProperty(Material::Density, Settings::instance()->setting("C/Density").toDouble());
    m_Materials[0].setProperty(Material::HeatCapacity, Settings::instance()->setting("C/HeatCapacity").toDouble());
    m_Materials[0].setProperty(Material::HeatConductivity, Settings::instance()->setting("C/HeatConductivity").toDouble());

    //****************************************************************************************************************************
    // Кремний (Расплав).
    //****************************************************************************************************************************
    m_Materials.append(Material(Material::Liquid));
    m_Materials[1].setProperty(Material::Diffusivity, 0);
    m_Materials[1].setProperty(Material::Density, Settings::instance()->setting("Si(Liquid)/Density").toDouble());
    m_Materials[1].setProperty(Material::MolarMass, Settings::instance()->setting("Si(Liquid)/MolarMass").toDouble());
    m_Materials[1].setProperty(Material::Viscosity, Settings::instance()->setting("Si(Liquid)/Viscosity").toDouble());
    m_Materials[1].setProperty(Material::HeatCapacity, Settings::instance()->setting("Si(Liquid)/HeatCapacity").toDouble());
    m_Materials[1].setProperty(Material::Expansiveness, Settings::instance()->setting("Si(Liquid)/Expansiveness").toDouble());
    m_Materials[1].setProperty(Material::HeatConductivity, Settings::instance()->setting("Si(Liquid)/HeatConductivity").toDouble());
    m_Materials[1].setProperty(Material::CrystallizationHeat, Settings::instance()->setting("Si(Liquid)/CrystallizationHeat").toDouble());

    //****************************************************************************************************************************
    // Кремний (Кристалл).
    //****************************************************************************************************************************
    m_Materials.append(Material(Material::Solid));
    m_Materials[2].setProperty(Material::Diffusivity, 0);
    m_Materials[2].setProperty(Material::Density, Settings::instance()->setting("Si(Solid)/Density").toDouble());
    m_Materials[2].setProperty(Material::HeatCapacity, Settings::instance()->setting("Si(Solid)/HeatCapacity").toDouble());
    m_Materials[2].setProperty(Material::HeatConductivity, Settings::instance()->setting("Si(Solid)/HeatConductivity").toDouble());

    //****************************************************************************************************************************
    // Аргон.
    //****************************************************************************************************************************
    m_Materials.append(Material(Material::Gaseous));
    m_Materials[3].setProperty(Material::Diffusivity, 0);
    m_Materials[3].setProperty(Material::Density, Settings::instance()->setting("Ar/Density").toDouble());
    m_Materials[3].setProperty(Material::HeatCapacity, Settings::instance()->setting("Ar/HeatCapacity").toDouble());
    m_Materials[3].setProperty(Material::HeatConductivity, Settings::instance()->setting("Ar/HeatConductivity").toDouble());

    //****************************************************************************************************************************
    // Фронт.
    //****************************************************************************************************************************
    m_Materials.append(Material(Material::Solid));
    m_Materials[4].setProperty(Material::Diffusivity, 0);
    m_Materials[4].setProperty(Material::Density, m_Materials[1].property(Material::Density));
    double capacity = m_Materials[1].property(Material::HeatCapacity) + m_Materials[1].property(Material::CrystallizationHeat)/(2*m_FrontHalfWidth);
    m_Materials[4].setProperty(Material::HeatCapacity, capacity);
    m_Materials[4].setProperty(Material::HeatConductivity, m_Materials[1].property(Material::HeatConductivity));

    //****************************************************************************************************************************
    // Примесь.
    //****************************************************************************************************************************
    m_Materials.append(Material(Material::Admixture));
    QString admixture = Settings::instance()->setting("Physics/Admixture").toString();
    m_Materials[5].setProperty(Material::MolarMass, Settings::instance()->setting(admixture + "/MolarMass").toDouble());
    m_Materials[5].setProperty(Material::Diffusivity, Settings::instance()->setting(admixture + "/Diffusivity").toDouble());
    m_Materials[5].setProperty(Material::SegregationFactor, Settings::instance()->setting(admixture + "/SegregationFactor").toDouble());

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    double density = m_Materials[1].property(Material::Density);
    double mass_c = m_Materials[1].property(Material::MolarMass);
    double mass_s = m_Materials[5].property(Material::MolarMass);
    double ppa = Settings::instance()->setting("Physics/Concentration").toDouble();

    double concentration = mass_s*density/(mass_c*1e6);

    m_C = linal::Vector(m_Net.numOfPoints(), concentration);
    m_COldTime = linal::Vector(m_Net.numOfPoints(), concentration);
}

//================================================================================================================================
//================================================================================================================================
QList<geometry::_2d::Point> CylinderTCt::identifyFront(double temperature)
{
    QList<geometry::_2d::Point> front_points;
    QList<geometry::_2d::Segment> front_segments = fem::triangles::isoline(m_Net, m_T, temperature);

    //****************************************************************************************************************************
    // Перевод списка сегментов в список точек.
    //****************************************************************************************************************************
    for(int i = 0; i < front_segments.size(); i++)
    {
        if(!front_points.contains(front_segments[i].begin()))
            front_points.append(front_segments[i].begin());
        if(!front_points.contains(front_segments[i].end()))
            front_points.append(front_segments[i].end());
    }

    //****************************************************************************************************************************
    // Удаление точек не из расплава.
    //****************************************************************************************************************************
    for(int i = 0; i < front_points.size(); ++i)
        if(front_points[i].x() > m_R[3] + DBL_EPSILON || front_points[i].y() + DBL_EPSILON < m_Z[3] || front_points[i].y() > m_Z[5] + DBL_EPSILON)
            front_points.removeAt(i--);

    //****************************************************************************************************************************
    // Сортировка точек по возрастанию координаты r.
    //****************************************************************************************************************************
    for(int i = front_points.size() - 1 ; i > 0 ; --i)
        for(int j = 0 ; j < i ; ++j)
            if(front_points[j].x() > front_points[j + 1].x())
                front_points.swap(j,j + 1);

    //****************************************************************************************************************************
    // Пристыковка фронта к стенкам тигля.
    //****************************************************************************************************************************
    if(!front_points.isEmpty())
    {
        if(qAbs(front_points.first().x()) < DBL_EPSILON)
            front_points.first().setX(0.0);
        if(qAbs(front_points.last().x() - m_R[3]) < DBL_EPSILON)
            front_points.last().setX(m_R[3]);
        if(qAbs(front_points.first().y() - m_Z[5]) < DBL_EPSILON)
            front_points.first().setY(m_Z[5]);
        if(qAbs(front_points.last().y() - m_Z[3]) < DBL_EPSILON)
            front_points.last().setY(m_Z[3]);
    }

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    for(int i = 0; i < front_points.size(); i++)
        if(qAbs(front_points[i].y() - m_Z[5]) < DBL_EPSILON)
        {
            for(int j = 0; j < i; j++)
                front_points.removeFirst();
            i = 0;
        }

    //****************************************************************************************************************************
    // Искуственный сдвиг фронта так, чтобы он не проходил точно через угол.
    //****************************************************************************************************************************
    if(!front_points.isEmpty())
    {
        if(qAbs(front_points.first().x() - m_R[4]) < DBL_EPSILON && qAbs(front_points.first().y() - m_Z[4]) < DBL_EPSILON)
            front_points.first().setX(m_R[4] + 10*DBL_EPSILON);

        if(qAbs(front_points.first().x() - m_R[3]) < DBL_EPSILON && qAbs(front_points.first().y() - m_Z[3]) < DBL_EPSILON)
            front_points.first().setY(m_Z[3] + 10*DBL_EPSILON);
    }

    return front_points;
}

//================================================================================================================================
//================================================================================================================================
QList<geometry::_2d::Point> CylinderTCt::identifyFrontExt(double temperature)
{
    QList<geometry::_2d::Point> front = identifyFront(temperature);

    if(front.isEmpty() || front.last().x() < 0.05)
    {
        double last_x;
        front.isEmpty() ? last_x = 0.0 : last_x = front.last().x();
        for(int i = 0; i <= (0.05 - last_x)/0.00125; ++i)
            front.append(geometry::_2d::Point(last_x + 0.00125*i, 0.005));
    }
    return front;
}

//================================================================================================================================
//================================================================================================================================
QVector<double> CylinderTCt::frontSpeed()
{
    QVector<double> speed(m_Front.size(), 0.0);

    for(int i = 1; i < speed.size() - 1; ++i)
    {
        geometry::_2d::Vector normal(m_Net.point(m_Front[i - 1]), m_Net.point(m_Front[i + 1]));
        normal = geometry::_2d::Vector::normal(normal);
        int triangle = m_NetSearcher->findTriangle(m_Net, m_Net.point(m_Front[i]).shifted(normal.x(), normal.y()));
        if(m_Materials[m_Net.material(triangle)].state() == Material::Liquid)
            normal = geometry::_2d::Vector(-normal.x(), -normal.y());

        for(int j = 0; j < m_OldFront.size() - 1; j++)
        {
            bool flag;
            double t,s;
            geometry::_2d::LineParam a(m_OldFront[j], m_OldFront[j + 1]);
            geometry::_2d::LineParam b(m_Net.point(m_Front[i]), normal);
            geometry::_2d::LineParam::findIntersection(a,b,&t,&s,&flag);
            if(!flag && t >= 0. && t <= 1.)
            {
                speed[i] = geometry::_2d::distance(m_Net.point(m_Front[i]), a.pos(t))/m_TimeStep;
                break;
            }
        }
    }
    speed[0] = speed[1];
    speed[speed.size() - 1] = speed[speed.size() - 2];

    return speed;
}

//================================================================================================================================
//================================================================================================================================
bool CylinderTCt::recreateNet()
{
    bool net_modified(false);
    fem::Triangulation triangulation(QStringList() << "T");

    //****************************************************************************************************************************
    //Формируем элементы базиса, которые не зависят от положения фронта
    //****************************************************************************************************************************
    int basis_points[13];
    basis_points[0] =  triangulation.addBasisPoint(geometry::_2d::Point(m_R[0],m_Z[0]));
    basis_points[1] =  triangulation.addBasisPoint(geometry::_2d::Point(m_R[1],m_Z[1]));
    basis_points[2] =  triangulation.addBasisPoint(geometry::_2d::Point(m_R[2],m_Z[2]));
    basis_points[3] =  triangulation.addBasisPoint(geometry::_2d::Point(m_R[3],m_Z[3]));
    basis_points[4] =  triangulation.addBasisPoint(geometry::_2d::Point(m_R[4],m_Z[4]));
    basis_points[5] =  triangulation.addBasisPoint(geometry::_2d::Point(m_R[5],m_Z[5]));
    basis_points[6] =  triangulation.addBasisPoint(geometry::_2d::Point(m_R[6],m_Z[6]));
    basis_points[7] =  triangulation.addBasisPoint(geometry::_2d::Point(m_R[7],m_Z[7]));
    basis_points[8] =  triangulation.addBasisPoint(geometry::_2d::Point(m_R[8],m_Z[8]));
    basis_points[9] =  triangulation.addBasisPoint(geometry::_2d::Point(m_R[9],m_Z[9]));
    basis_points[10] = triangulation.addBasisPoint(geometry::_2d::Point(m_R[10],m_Z[10]));
    basis_points[11] = triangulation.addBasisPoint(geometry::_2d::Point(m_R[11],m_Z[11]));
    basis_points[12] = triangulation.addBasisPoint(geometry::_2d::Point(m_R[12],m_Z[12]));

    int border_segments[19];
    border_segments[0] =  triangulation.addBorderSegment(basis_points[0], basis_points[1], m_N[0]);
    border_segments[1] =  triangulation.addBorderSegment(basis_points[0], basis_points[2], m_N[1]);
    border_segments[2] =  triangulation.addBorderSegment(basis_points[1], basis_points[3], m_N[2]);
    border_segments[6] =  triangulation.addBorderSegment(basis_points[1], basis_points[6], m_N[6]);
    border_segments[8] =  triangulation.addBorderSegment(basis_points[5], basis_points[6], m_N[8]);
    border_segments[9] =  triangulation.addBorderSegment(basis_points[4], basis_points[7], m_N[9]);
    border_segments[10] = triangulation.addBorderSegment(basis_points[5], basis_points[8], m_N[10]);
    border_segments[11] = triangulation.addBorderSegment(basis_points[6], basis_points[10],m_N[11]);
    border_segments[12] = triangulation.addBorderSegment(basis_points[7], basis_points[8], m_N[12]);
    border_segments[13] = triangulation.addBorderSegment(basis_points[8], basis_points[10],m_N[13]);
    border_segments[14] = triangulation.addBorderSegment(basis_points[7], basis_points[11],m_N[14]);
    border_segments[15] = triangulation.addBorderSegment(basis_points[7], basis_points[9], m_N[15]);
    border_segments[16] = triangulation.addBorderSegment(basis_points[9], basis_points[10],m_N[16]);
    border_segments[17] = triangulation.addBorderSegment(basis_points[9], basis_points[12],m_N[17]);
    border_segments[18] = triangulation.addBorderSegment(basis_points[11],basis_points[12],m_N[18]);

    triangulation.addBorderSegmentBCInfo("T", border_segments[0], 1, 0);
    triangulation.addBorderSegmentBCInfo("T", border_segments[6], 1, 1);
    triangulation.addBorderSegmentBCInfo("T", border_segments[11], 1, 1);
    triangulation.addBorderSegmentBCInfo("T", border_segments[18], 1, m_TRod != 0 ? 2 : 1);

    triangulation.addQuadrangleSubdomain(border_segments[13], border_segments[11], border_segments[8],  border_segments[10],0);
    triangulation.addQuadrangleSubdomain(border_segments[16], border_segments[13], border_segments[12], border_segments[15],0);
    triangulation.addQuadrangleSubdomain(border_segments[18], border_segments[17], border_segments[15], border_segments[14],0);

    //****************************************************************************************************************************
    //****************************************************************************************************************************
    //Изотермы фронта.
    QList<geometry::_2d::Point> front_plus_points = identifyFront(m_TCrystal + m_FrontHalfWidth);
    QList<geometry::_2d::Point> front_minus_points = identifyFront(m_TCrystal - m_FrontHalfWidth);

    QList<double> front_plus_r, front_plus_z;
    QList<double> front_minus_r, front_minus_z;

    for(int i = 0; i < front_plus_points.size(); ++i)
    {
        front_plus_r.append(front_plus_points[i].x());
        front_plus_z.append(front_plus_points[i].y());
    }

    for(int i = 0; i < front_minus_points.size(); ++i)
    {
        front_minus_r.append(front_minus_points[i].x());
        front_minus_z.append(front_minus_points[i].y());
    }

    //HACK: Кристалл очень маленький. Считаем, что его нет.
    if((front_minus_points.empty() || front_minus_points.first().y() <= m_Z[2] + m_MinCrystalHeight) && m_State == SolverTCt::NotStarted)
    {
        front_minus_points.clear();
        front_plus_points.clear();
    }

    //****************************************************************************************************************************
    //Нету фронта кристаллизации
    //****************************************************************************************************************************
    if(front_plus_points.empty() && front_minus_points.empty())
    {
        if(m_State == SolverTCt::NotStarted || m_State == SolverTCt::Finished)
            net_modified = false;
        else
        {
            border_segments[3] =  triangulation.addBorderSegment(basis_points[2], basis_points[3], m_N[3]);
            border_segments[4] =  triangulation.addBorderSegment(basis_points[2], basis_points[4], m_N[4]);
            border_segments[5] =  triangulation.addBorderSegment(basis_points[3], basis_points[5], m_N[5]);
            border_segments[7] =  triangulation.addBorderSegment(basis_points[4], basis_points[5], m_N[7]);

            triangulation.addBorderSegmentBCInfo("T", border_segments[3], 4);
            triangulation.addBorderSegmentBCInfo("T", border_segments[4], 2, 0);
            triangulation.addBorderSegmentBCInfo("T", border_segments[5], 4);
            triangulation.addBorderSegmentBCInfo("T", border_segments[7], 4);

            triangulation.addQuadrangleSubdomain(border_segments[3],  border_segments[2],  border_segments[0],  border_segments[1], 0);
            triangulation.addQuadrangleSubdomain(border_segments[8],  border_segments[6],  border_segments[2],  border_segments[5], 0);
            triangulation.addQuadrangleSubdomain(border_segments[12], border_segments[10], border_segments[7],  border_segments[9], 3);

            if(m_State == SolverTCt::JustStarted)
            {
                triangulation.addQuadrangleSubdomain(border_segments[7],  border_segments[5],  border_segments[3],  border_segments[4], 1);
                m_State = SolverTCt::NotStarted;
            }
            else if(m_State == SolverTCt::InProgress)
            {

                if(m_T[basis_points[4]] < m_TCrystal)
                {
                    triangulation.addQuadrangleSubdomain(border_segments[7],  border_segments[5],  border_segments[3],  border_segments[4], 2);
                    m_State = SolverTCt::Finished;
                }
                else
                {
                    triangulation.addQuadrangleSubdomain(border_segments[7],  border_segments[5],  border_segments[3],  border_segments[4], 1);
                    m_State = SolverTCt::NotStarted;
                }
            }
            net_modified = true;
        }
    }

    //****************************************************************************************************************************
    //Есть обе границы фронта
    //****************************************************************************************************************************
    else if(!front_minus_points.empty() && !front_plus_points.empty())
    {
        int added_points[4];
        added_points[0] = triangulation.addBasisPoint(front_minus_points.first());
        added_points[1] = triangulation.addBasisPoint(front_minus_points.last());
        added_points[2] = triangulation.addBasisPoint(front_plus_points.first());
        added_points[3] = triangulation.addBasisPoint(front_plus_points.last());

        double r_plus_first(front_plus_points.first().x()), r_plus_last(front_plus_points.last().x());
        double z_plus_first(front_plus_points.first().y()), z_plus_last(front_plus_points.last().y());
        double r_minus_first(front_minus_points.first().x()), r_minus_last(front_minus_points.last().x());
        double z_minus_first(front_minus_points.first().y()), z_minus_last(front_minus_points.last().y());

        QList<int> front_plus_indexes;
        QList<int> front_minus_indexes;
        // Распределяем точки на фронте равномерно по z.
        if(r_minus_first > 0)
        {
            QList< QVector<double> > front_plus_interpolated = approximation::Interpolation::interpolate(front_plus_z, front_plus_r, m_N[5]*(m_Z[5] - z_plus_last)/m_L[5], approximation::Interpolation::Linear);
            QList< QVector<double> > front_minus_interpolated = approximation::Interpolation::interpolate(front_minus_z, front_minus_r, m_N[5]*(m_Z[5] - z_minus_last)/m_L[5], approximation::Interpolation::Linear);

            front_plus_indexes.append(added_points[2]);
            for(int i = 0; i < front_plus_interpolated.size(); i++)
                front_plus_indexes.append(triangulation.addBasisPoint(geometry::_2d::Point(front_plus_interpolated[i][1],front_plus_interpolated[i][0])));
            front_plus_indexes.append(added_points[3]);

            front_minus_indexes.append(added_points[0]);
            for(int i = 0; i < front_minus_interpolated.size(); i++)
                front_minus_indexes.append(triangulation.addBasisPoint(geometry::_2d::Point(front_minus_interpolated[i][1],front_minus_interpolated[i][0])));
            front_minus_indexes.append(added_points[1]);
        }
        // или равномерно по r.
        else
        {
            QList< QVector<double> > front_plus_interpolated = approximation::Interpolation::interpolate(front_plus_r, front_plus_z, m_N[3]*r_plus_last/m_L[3], approximation::Interpolation::Linear);
            QList< QVector<double> > front_minus_interpolated = approximation::Interpolation::interpolate(front_minus_r, front_minus_z, m_N[3]*r_minus_last/m_L[3], approximation::Interpolation::Linear);

            front_plus_indexes.append(added_points[2]);
            for(int i = 0; i < front_plus_interpolated.size(); i++)
                front_plus_indexes.append(triangulation.addBasisPoint(geometry::_2d::Point(front_plus_interpolated[i][0],front_plus_interpolated[i][1])));
            front_plus_indexes.append(added_points[3]);

            front_minus_indexes.append(added_points[0]);
            for(int i = 0; i < front_minus_interpolated.size(); i++)
                front_minus_indexes.append(triangulation.addBasisPoint(geometry::_2d::Point(front_minus_interpolated[i][0],front_minus_interpolated[i][1])));
            front_minus_indexes.append(added_points[1]);
        }

        m_Front = front_plus_indexes;

        //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        // Левая граница фронта находится на оси симметрии, правая - на дне тигля.
        //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        if(r_plus_last < m_R[3])
        {
            int added_segments[12];
            added_segments[0] = triangulation.addBorderSegment(basis_points[2],added_points[1],front_minus_indexes.size() - 2);
            added_segments[1] = triangulation.addBorderSegment(added_points[3],basis_points[3],(int)(m_N[3]*(m_R[3] - r_plus_last)/m_L[3]));
            added_segments[2] = triangulation.addBorderSegment(basis_points[2],added_points[0],(int)(m_N[4]*(z_minus_first - m_Z[2])/m_L[4]));
            added_segments[3] = triangulation.addBorderSegment(added_points[2],basis_points[4],(int)(m_N[4]*(m_Z[4] - z_plus_first)/m_L[4]));
            added_segments[4] = triangulation.addBorderSegment(fem::BorderSegment(front_minus_indexes));
            added_segments[5] = triangulation.addBorderSegment(fem::BorderSegment(front_plus_indexes));
            added_segments[6] = triangulation.addBorderSegment(added_points[1],added_points[3],front_plus_indexes.size() - front_minus_indexes.size());
            added_segments[7] = triangulation.addBorderSegment(added_points[0],added_points[2],0);
            added_segments[8] = triangulation.concatinateBorderSegments(added_segments[0],added_segments[6]);
            added_segments[9] = triangulation.concatinateBorderSegments(added_segments[8],added_segments[1]);
            added_segments[10] = triangulation.concatinateBorderSegments(added_segments[5],added_segments[1]);
            added_segments[11] = triangulation.concatinateBorderSegments(added_segments[4],added_segments[6]);

            border_segments[5] = triangulation.addBorderSegment(basis_points[3],basis_points[5],m_N[5]);
            border_segments[7] = triangulation.addBorderSegment(basis_points[4],basis_points[5],triangulation.pointsCount(added_segments[10]) - 2);

            triangulation.addQuadrangleSubdomain(border_segments[8], border_segments[6], border_segments[2], border_segments[5], 0);
            triangulation.addQuadrangleSubdomain(border_segments[12], border_segments[10], border_segments[7], border_segments[9], 3);
            triangulation.addQuadrangleSubdomain(added_segments[9] ,border_segments[2], border_segments[0], border_segments[1], 0);
            triangulation.addQuadrangleSubdomain(border_segments[7], border_segments[5], added_segments[10],added_segments[3], 1);
            triangulation.addTriangleSubdomain(added_segments[0], added_segments[4], added_segments[2], 2);
            triangulation.addTriangleSubdomain(added_segments[5],added_segments[11],added_segments[7],4);
        }
        //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        // Левая граница находится на оси симметрии, правая граница фронта "размазана" по разные стороны угла тигля (у дна)
        //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        else if(z_plus_last > m_Z[3] && r_minus_last < m_R[3])
        {
            border_segments[7] = triangulation.addBorderSegment(basis_points[4],basis_points[5],front_plus_indexes.size() - 2);

            int added_segments[12];
            added_segments[0] = triangulation.addBorderSegment(basis_points[2],added_points[1],front_minus_indexes.size() - 1);
            added_segments[1] = triangulation.addBorderSegment(added_points[1],basis_points[3],m_N[3] - front_minus_indexes.size());
            added_segments[2] = triangulation.addBorderSegment(basis_points[3],added_points[3],0);
            added_segments[3] = triangulation.addBorderSegment(added_points[3],basis_points[5],m_N[5] - 1);
            added_segments[4] = triangulation.addBorderSegment(fem::BorderSegment(front_minus_indexes));
            added_segments[5] = triangulation.addBorderSegment(fem::BorderSegment(front_plus_indexes));
            added_segments[6] = triangulation.addBorderSegment(added_points[0],added_points[2],0);
            added_segments[7] = triangulation.addBorderSegment(basis_points[2],added_points[0],(int)(m_N[4]*(z_minus_first - m_Z[2])/m_L[4]));
            added_segments[8] = triangulation.addBorderSegment(added_points[2],basis_points[4],(int)(m_N[4]*(m_Z[4] - z_plus_first)/m_L[4]));
            added_segments[9] = triangulation.concatinateBorderSegments(added_segments[2],added_segments[3]);
            added_segments[10] = triangulation.concatinateBorderSegments(added_segments[0],added_segments[1]);
            added_segments[11] = triangulation.concatinateBorderSegments(added_segments[4],added_segments[1]);

            triangulation.addQuadrangleSubdomain(border_segments[8], border_segments[6], border_segments[2], added_segments[9], 0);
            triangulation.addQuadrangleSubdomain(border_segments[12], border_segments[10], border_segments[7], border_segments[9], 3);
            triangulation.addQuadrangleSubdomain(added_segments[10],border_segments[2], border_segments[0],border_segments[1], 0);
            triangulation.addQuadrangleSubdomain(border_segments[7], added_segments[3], added_segments[5],added_segments[8], 1);
            triangulation.addTriangleSubdomain(added_segments[0], added_segments[4], added_segments[7], 2);
            triangulation.addQuadrangleSubdomain(added_segments[5],added_segments[2],added_segments[11],added_segments[6],4);
        }
        //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        // Левая граница фронта на оси симметрии, правая - на боковой стенке тигля.
        //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        else if(z_plus_first < m_Z[4])
        {
            border_segments[7] = triangulation.addBorderSegment(basis_points[4],basis_points[5],front_plus_indexes.size() - 2);
            border_segments[3] = triangulation.addBorderSegment(basis_points[2],basis_points[3],front_minus_indexes.size() - 2);

            int added_segments[10];
            added_segments[0] = triangulation.addBorderSegment(basis_points[3],added_points[1],(int)(m_N[5]*(z_minus_last - m_Z[3])/m_L[5]));
            added_segments[1] = triangulation.addBorderSegment(added_points[3],basis_points[5],(int)(m_N[5]*(m_Z[5] - z_plus_last)/m_L[5]));
            added_segments[2] = triangulation.addBorderSegment(basis_points[2],added_points[0],(int)(m_N[4]*(z_minus_first - m_Z[2])/m_L[4]));
            added_segments[3] = triangulation.addBorderSegment(added_points[2],basis_points[4],(int)(m_N[4]*(m_Z[4] - z_plus_first)/m_L[4]));
            added_segments[4] = triangulation.addBorderSegment(fem::BorderSegment(front_minus_indexes));
            added_segments[5] = triangulation.addBorderSegment(fem::BorderSegment(front_plus_indexes));
            added_segments[6] = triangulation.addBorderSegment(added_points[0],added_points[2],0);
            added_segments[7] = triangulation.addBorderSegment(added_points[1],added_points[3],0);
            added_segments[8] = triangulation.concatinateBorderSegments(added_segments[0],added_segments[7]);
            added_segments[9] = triangulation.concatinateBorderSegments(added_segments[8],added_segments[1]);

            triangulation.addQuadrangleSubdomain(border_segments[3], border_segments[2], border_segments[0], border_segments[1], 0);
            triangulation.addQuadrangleSubdomain(border_segments[12], border_segments[10], border_segments[7], border_segments[9], 3);
            triangulation.addQuadrangleSubdomain(border_segments[8], border_segments[6], border_segments[2], added_segments[9], 0);
            triangulation.addQuadrangleSubdomain(border_segments[7], added_segments[1], added_segments[5], added_segments[3], 1);
            triangulation.addQuadrangleSubdomain(added_segments[4], added_segments[0], border_segments[3], added_segments[2], 2);
            triangulation.addQuadrangleSubdomain(added_segments[5], added_segments[7], added_segments[4], added_segments[6], 4);
        }
        //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        // Левая граница фронта "размазана" по разные стороны угла (у свободной поверхности), правая - на боковой стенке.
        //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        else if(r_plus_first > m_R[4] && z_minus_first < m_Z[4])
        {
            border_segments[3] = triangulation.addBorderSegment(basis_points[2],basis_points[3], front_minus_indexes.size() - 2);

            int added_segments[13];
            added_segments[0] = triangulation.addBorderSegment(basis_points[2],added_points[0],(int)(m_N[4]*(z_minus_first - m_Z[2])/m_L[4]));
            added_segments[1] = triangulation.addBorderSegment(added_points[0],basis_points[4],0);
            added_segments[2] = triangulation.addBorderSegment(basis_points[4],added_points[2],0);
            added_segments[3] = triangulation.addBorderSegment(added_points[2],basis_points[5],(int)(m_N[7]*(m_R[5] - r_plus_first)/m_L[5]));
            added_segments[4] = triangulation.addBorderSegment(basis_points[3],added_points[1],(int)(m_N[5]*(z_minus_last - m_Z[3])/m_L[5]));
            added_segments[5] = triangulation.addBorderSegment(added_points[3],basis_points[5],front_plus_indexes.size() - 2);
            added_segments[6] = triangulation.addBorderSegment(fem::BorderSegment(front_minus_indexes));
            added_segments[7] = triangulation.addBorderSegment(fem::BorderSegment(front_plus_indexes));
            added_segments[8] = triangulation.addBorderSegment(added_points[1],added_points[3],0);
            added_segments[9] = triangulation.concatinateBorderSegments(added_segments[2],added_segments[3]);
            added_segments[10] = triangulation.concatinateBorderSegments(added_segments[4],added_segments[8]);
            added_segments[11] = triangulation.concatinateBorderSegments(added_segments[10],added_segments[5]);
            added_segments[12] = triangulation.concatinateBorderSegments(added_segments[2], added_segments[7]);

            triangulation.addQuadrangleSubdomain(border_segments[3], border_segments[2], border_segments[0], border_segments[1], 0);
            triangulation.addQuadrangleSubdomain(added_segments[6], added_segments[4], border_segments[3], added_segments[0], 2);
            triangulation.addQuadrangleSubdomain(border_segments[12], border_segments[10], added_segments[9], border_segments[9], 3);
            triangulation.addQuadrangleSubdomain(border_segments[8], border_segments[6], border_segments[2], added_segments[11], 0);
            triangulation.addTriangleSubdomain(added_segments[3], added_segments[7], added_segments[5], 1);
            triangulation.addQuadrangleSubdomain(added_segments[12], added_segments[8], added_segments[6], added_segments[1], 4);
        }
        //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        // Левая граница фронта лежит на свободной поверхности, правая - на боковой стенке тигля.
        //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        else if(r_plus_first < m_R[5])
        {
            border_segments[3] = triangulation.addBorderSegment(basis_points[2],basis_points[3],front_minus_indexes.size() + (int)(m_N[7]*r_minus_first/m_L[7]));
            border_segments[4] = triangulation.addBorderSegment(basis_points[2],basis_points[4],m_N[5] - 1);

            fem::BorderSegment front_minus_segment(front_minus_indexes);

            int added_segments[14];
            added_segments[0] = triangulation.addBorderSegment(basis_points[4],added_points[0],(int)(m_N[7]*r_minus_first/m_L[7]));
            added_segments[1] = triangulation.addBorderSegment(added_points[2],basis_points[5],(int)(m_N[7]*(m_R[5] - r_plus_first)/m_L[7]));
            added_segments[2] = triangulation.addBorderSegment(basis_points[3],added_points[1],(int)(m_N[5]*(z_minus_last - m_Z[3])/m_L[5]));
            added_segments[3] = triangulation.addBorderSegment(added_points[3],basis_points[5],front_plus_indexes.size() - 2);
            added_segments[4] = triangulation.addBorderSegment(front_minus_segment);
            front_minus_segment.reverse();
            added_segments[13] = triangulation.addBorderSegment(front_minus_segment);
            added_segments[5] = triangulation.addBorderSegment(fem::BorderSegment(front_plus_indexes));
            added_segments[6] = triangulation.addBorderSegment(added_points[0],added_points[2],0);
            added_segments[7] = triangulation.addBorderSegment(added_points[1],added_points[3],0);
            added_segments[8] = triangulation.concatinateBorderSegments(added_segments[2],added_segments[13]);
            added_segments[9] = triangulation.concatinateBorderSegments(added_segments[0],added_segments[6]);
            added_segments[10] = triangulation.concatinateBorderSegments(added_segments[9],added_segments[1]);
            added_segments[11] = triangulation.concatinateBorderSegments(added_segments[2], added_segments[7]);
            added_segments[12] = triangulation.concatinateBorderSegments(added_segments[11],added_segments[3]);

            triangulation.addQuadrangleSubdomain(border_segments[3], border_segments[2], border_segments[0], border_segments[1], 0);
            triangulation.addQuadrangleSubdomain(added_segments[0], added_segments[8], border_segments[3], border_segments[4], 2);
            triangulation.addQuadrangleSubdomain(border_segments[12], border_segments[10], added_segments[10], border_segments[9], 3);
            triangulation.addQuadrangleSubdomain(border_segments[8], border_segments[6], border_segments[2],added_segments[12], 0);
            triangulation.addTriangleSubdomain(added_segments[3], added_segments[5], added_segments[1], 1);
            triangulation.addQuadrangleSubdomain(added_segments[5],added_segments[7],added_segments[4],added_segments[6], 4);
        }
        m_State = SolverTCt::InProgress;
        net_modified = true;
    }

    //****************************************************************************************************************************
    //Собственно перестройка сетки и все операции которые из нее следуют.
    //****************************************************************************************************************************
    if(net_modified)
    {
        fem::TriangleNet new_net = triangulation.triangulate();
        triangulation.getDirichletBoundaryInfo("T", m_TopTBC, m_TypeTBC);

        //------------------------------------------------------------------------------------------------------------------------
        //------------------------------------------------------------------------------------------------------------------------
        static linal::Vector new_t;
        static linal::Vector new_t_old_time;

        static linal::Vector new_c;
        static linal::Vector new_c_old_time;

        new_t.resize(new_net.numOfPoints());
        new_t_old_time.resize(new_net.numOfPoints());

        new_c.resize(new_net.numOfPoints());
        new_c_old_time.resize(new_net.numOfPoints());

        for(int i = 0; i < new_net.numOfPoints(); i++)
        {
            int triangle = m_NetSearcher->findTriangle(m_Net, new_net.point(i));

            new_t[i] = fem::triangles::functionValue(triangle, m_Net, m_T, new_net.point(i));
            new_t_old_time[i] = fem::triangles::functionValue(triangle, m_Net, m_TOldTime, new_net.point(i));

            new_c[i] = fem::triangles::functionValue(triangle, m_Net, m_C, new_net.point(i));
            new_c_old_time[i] = fem::triangles::functionValue(triangle, m_Net, m_COldTime, new_net.point(i));
        }

        m_T = new_t;
        m_TOldTime = new_t_old_time;

        m_C = new_c;
        m_COldTime = new_c_old_time;

        //------------------------------------------------------------------------------------------------------------------------
        //------------------------------------------------------------------------------------------------------------------------
        m_Net = new_net;
        m_NetSearcher->analyzeNet(m_Net);

        //------------------------------------------------------------------------------------------------------------------------
        //------------------------------------------------------------------------------------------------------------------------
        m_F.resize(m_Net.numOfPoints());

        QVector<int> igl, jgl, igu, ngu;
        fem::triangles::formMatrixProfile(m_Net, igl, jgl, igu, ngu);
        m_A = linal::CSRMatrix(m_Net.numOfPoints(), igl, jgl, igu, ngu);

        foreach(linal::SLAESolver* solver, m_Solvers)
            solver->adaptForSLAE(m_A.dimension());

        //------------------------------------------------------------------------------------------------------------------------
        //------------------------------------------------------------------------------------------------------------------------
        return true;
    }
    else
        return false;
}
